Analyse Individual Dynos Captured with Stereo-seq¶

In [1]:
import stereo as st
import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib as plt
import sys
import os
import matplotlib.pyplot as plt
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Information on chip:

  • horizontal dynos: Control + chir
  • vertial dynos: DYMK + chir

Load in chip at bin20 then section different dynos for analysis¶

In [2]:
## Load into stereoexpdata object
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
                                 spatial_key='spatial')
dyno_dm
Out[2]:
StereoExpData object with n_cells X n_genes = 13099 X 42359
bin_type: None
offset_x = None
offset_y = None
cells: ['cell_name', 'total_counts', 'pct_counts_mt', 'n_genes_by_counts']
genes: ['gene_name']
cells_matrix = ['spatial']
Layers with keys: 
tl.result: []
In [3]:
ins = dyno_dm.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
ins.show()
Out[3]:
In [4]:
dyno_dm_1 = dyno_dm.tl.filter_coordinates(min_x=9600,max_x=12600,min_y=12600,max_y=14600)
ins = dyno_dm_1.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)

dyno_dm_1.tl.filter_cells(
        min_gene=200,
        min_n_genes_by_counts=3,
        pct_counts_mt=5,
        inplace=True
        )
dyno_dm_1 = st.io.stereo_to_anndata(dyno_dm_1,flavor='seurat',
                                output='dyno_dm_1.h5ad')
ins.show()
[2025-04-14 14:47:10][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates...
[2025-04-14 14:47:10][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0054s.
[2025-04-14 14:47:13][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells...
[2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0404s.
[2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident'].
[2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] .
[2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] .
[2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info.
[2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata.
[2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_1.h5ad
Out[4]:
In [5]:
## Load into stereoexpdata object
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
                                 spatial_key='spatial')

dyno_dm_2 = dyno_dm.tl.filter_coordinates(min_x=10000,max_x=11500,min_y=5700,max_y=9300)
ins = dyno_dm_2.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)

dyno_dm_2.tl.filter_cells(
        min_gene=200,
        min_n_genes_by_counts=3,
        pct_counts_mt=5,
        inplace=True
        )
dyno_dm_2 = st.io.stereo_to_anndata(dyno_dm_2,flavor='seurat',
                                output='dyno_dm_2.h5ad')
ins.show()
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates...
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0013s.
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells...
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0373s.
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident'].
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] .
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] .
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info.
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata.
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_2.h5ad
Out[5]:
In [6]:
## Load into stereoexpdata object
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
                                 spatial_key='spatial')

dyno_dm_3 =dyno_dm.tl.filter_coordinates(min_x=13600,max_x=16600,min_y=16600,max_y=19600)
ins = dyno_dm_3.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
dyno_dm_3.tl.filter_cells(
        min_gene=200,
        min_n_genes_by_counts=3,
        pct_counts_mt=5,
        inplace=True
        )

dyno_dm_3 = st.io.stereo_to_anndata(dyno_dm_3,flavor='seurat',
                                output='dyno_dm_3.h5ad')
ins.show()
[2025-04-14 14:47:16][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates...
[2025-04-14 14:47:16][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0014s.
[2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells...
[2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0761s.
[2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident'].
[2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] .
[2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] .
[2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info.
[2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata.
[2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_3.h5ad
Out[6]:
In [7]:
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
                                 spatial_key='spatial')

dyno_dm_4 = dyno_dm.tl.filter_coordinates(min_x=9600,max_x=12600,min_y=16000,max_y=18600)
dyno_dm_4
dyno_dm_4.tl.filter_cells(
        min_gene=200,
        min_n_genes_by_counts=3,
        pct_counts_mt=5,
        inplace=True
        )
ins = dyno_dm_4.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)

dyno_dm_4 = st.io.stereo_to_anndata(dyno_dm_4,flavor='seurat',
                                output='dyno_dm_4.h5ad')
ins.show()
[2025-04-14 14:47:18][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates...
[2025-04-14 14:47:18][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0015s.
[2025-04-14 14:47:18][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells...
[2025-04-14 14:47:18][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0763s.
[2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident'].
[2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] .
[2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] .
[2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info.
[2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata.
[2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_4.h5ad
Out[7]:
In [8]:
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
                                 spatial_key='spatial')
dyno_dm_5= dyno_dm.tl.filter_coordinates(min_x=11700,max_x=13200,min_y=3200,max_y=10150)
dyno_dm_5.tl.filter_cells(
        min_gene=200,
        min_n_genes_by_counts=3,
        pct_counts_mt=5,
        inplace=True
        )
ins = dyno_dm_5.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)

dyno_dm_5 = st.io.stereo_to_anndata(dyno_dm_5,flavor='seurat',
                                output='dyno_dm_5.h5ad')
ins.show()
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates...
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0013s.
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells...
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0363s.
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident'].
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] .
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] .
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info.
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata.
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_5.h5ad
Out[8]:
In [9]:
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
                                 spatial_key='spatial')
dyno_dm_6= dyno_dm.tl.filter_coordinates(min_x=13000,max_x=15000,min_y=3400,max_y=10500)
ins = dyno_dm_6.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)

dyno_dm_6.tl.filter_cells(
        min_gene=200,
        min_n_genes_by_counts=3,
        pct_counts_mt=5,
        inplace=True
        )
dyno_dm_6 = st.io.stereo_to_anndata(dyno_dm_6,flavor='seurat',
                                output='dyno_dm_6.h5ad')
ins.show()
[2025-04-14 14:47:21][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates...
[2025-04-14 14:47:21][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0013s.
[2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells...
[2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0804s.
[2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident'].
[2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] .
[2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] .
[2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info.
[2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata.
[2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_6.h5ad
Out[9]:
In [10]:
import anndata as ad
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
from shapely.geometry import Polygon, Point
#import decoupler as dc
import scanpy as sc
from sklearn.neighbors import NearestNeighbors
from scipy import stats
import squidpy as sq
#from pydeseq2.dds import DeseqDataSet, DefaultInference
#from pydeseq2.ds import DeseqStats
In [11]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA

adata1=ad.read('dyno_dm_1.h5ad')

# Extract spatial coordinates
points = adata1.obsm['spatial']

# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)

# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
    point = Point(row['x'], row['y'])
    return point.distance(polygon.boundary)

# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
    compute_distance_to_boundary, polygon=hull_polygon, axis=1
)

# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()

# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']

# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0]  # The first principal component represents the longer axis

# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis

# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj

# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3

# Classification function
def classify_region(projection, threshold_1, threshold_2):
    if projection <= threshold_1 or projection >= threshold_2:
        return 'pole'
    else:
        return 'center'

# Apply classification on filtered points
adata_filtered.obs['region'] = [
    classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')

# Visualize results
sq.pl.spatial_scatter(
    adata_filtered,
    size=10,
    library_id="spatial",
    figsize=(5, 5),
    frameon=False,
    shape=None,
    color=["region"]
)
adata_filtered.write('dyno1.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  _cax = scatter(
No description has been provided for this image
In [8]:
from matplotlib.legend import Legend 
In [12]:
adata_filtered=ad.read('dyno1.h5ad')
    
sq.pl.spatial_scatter(
    adata_filtered,
    size=10,
    library_id="spatial",
    figsize=(5, 5),
    frameon=False,
    shape=None,
    color=["region"],title=''
)
fig = plt.gcf()

# Remove the legend
for ax in fig.get_axes():
    for obj in ax.get_children():
        if isinstance(obj, Legend):
            obj.remove()

# **Manually Find and Remove the Colorbar**
for ax in fig.axes:
    if "colorbar" in str(ax):  # Check if this is the colorbar axis
        fig.delaxes(ax)

# **Second Attempt: Remove Extra Axes (Colorbars are usually last axes added)**
while len(fig.axes) > 1:  # Assuming 1 main plot, remove extras
    fig.delaxes(fig.axes[-1])  # Delete last added axis (likely colorbar)
plt.savefig('/group/tran3/mnucera/STOmics_crossorganoid_project/STOmics_notebooks/figures_stomics/fig4A.svg',dpi=300)
plt.savefig('/group/tran3/mnucera/STOmics_crossorganoid_project/STOmics_notebooks/figures_stomics/fig4A.svg',dpi=300)
plt.show()
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  _cax = scatter(
No description has been provided for this image
In [13]:
adata_filtered=ad.read('dyno5.h5ad')
    
sq.pl.spatial_scatter(
    adata_filtered,
    size=10,
    library_id="spatial",
    figsize=(5, 5),
    frameon=False,
    shape=None,
    color=["region"],title=''
)
fig = plt.gcf()

# Remove the legend
for ax in fig.get_axes():
    for obj in ax.get_children():
        if isinstance(obj, Legend):
            obj.remove()

# **Manually Find and Remove the Colorbar**
for ax in fig.axes:
    if "colorbar" in str(ax):  # Check if this is the colorbar axis
        fig.delaxes(ax)

# **Second Attempt: Remove Extra Axes (Colorbars are usually last axes added)**
while len(fig.axes) > 1:  # Assuming 1 main plot, remove extras
    fig.delaxes(fig.axes[-1])  # Delete last added axis (likely colorbar)
plt.savefig('/group/tran3/mnucera/STOmics_crossorganoid_project/STOmics_notebooks/figures_stomics/fig4B.svg',dpi=300)

plt.show()
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  _cax = scatter(
No description has been provided for this image
In [3]:
def create_custom_legend(labels, colors, filename="legend.svg", dpi=300):
    fig, ax = plt.subplots()
    ax.axis("off")  # Hide axes

    # Create legend handles
    handles = [plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=color, markersize=10) 
               for color in colors]

    # Add legend
    legend = ax.legend(handles, labels, loc="center", frameon=False)

    # Save figure with tight bounding box to remove excess whitespace
    fig.savefig(filename, dpi=dpi, bbox_inches="tight", transparent=True)

    plt.show()

# Example usage
labels = ["Center","Poles"]
colors = ['blue','orange']

create_custom_legend(labels, colors, filename="/group/tran3/mnucera/STOmics_crossorganoid_project/STOmics_notebooks/figures_stomics/dyno_legend.svg")
No description has been provided for this image
In [16]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA

adata1=ad.read('dyno_dm_2.h5ad')

# Extract spatial coordinates
points = adata1.obsm['spatial']

# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)

# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
    point = Point(row['x'], row['y'])
    return point.distance(polygon.boundary)

# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
    compute_distance_to_boundary, polygon=hull_polygon, axis=1
)

# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()

# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']

# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0]  # The first principal component represents the longer axis

# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis

# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj

# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3

# Classification function
def classify_region(projection, threshold_1, threshold_2):
    if projection <= threshold_1 or projection >= threshold_2:
        return 'pole'
    else:
        return 'center'

# Apply classification on filtered points
adata_filtered.obs['region'] = [
    classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')

# Visualize results
sq.pl.spatial_scatter(
    adata_filtered,
    size=10,
    library_id="spatial",
    figsize=(5, 5),
    frameon=False,
    shape=None,
    color=["region"]
)
adata_filtered.write('dyno2.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  _cax = scatter(
No description has been provided for this image
In [17]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA

adata1=ad.read('dyno_dm_3.h5ad')

# Extract spatial coordinates
points = adata1.obsm['spatial']

# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)

# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
    point = Point(row['x'], row['y'])
    return point.distance(polygon.boundary)

# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
    compute_distance_to_boundary, polygon=hull_polygon, axis=1
)

# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()

# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']

# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0]  # The first principal component represents the longer axis

# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis

# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj

# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3

# Classification function
def classify_region(projection, threshold_1, threshold_2):
    if projection <= threshold_1 or projection >= threshold_2:
        return 'pole'
    else:
        return 'center'

# Apply classification on filtered points
adata_filtered.obs['region'] = [
    classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')

# Visualize results
sq.pl.spatial_scatter(
    adata_filtered,
    size=10,
    library_id="spatial",
    figsize=(5, 5),
    frameon=False,
    shape=None,
    color=["region"]
)
adata_filtered.write('dyno3.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  _cax = scatter(
No description has been provided for this image
In [18]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA

adata1=ad.read('dyno_dm_4.h5ad')

# Extract spatial coordinates
points = adata1.obsm['spatial']

# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)

# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
    point = Point(row['x'], row['y'])
    return point.distance(polygon.boundary)

# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
    compute_distance_to_boundary, polygon=hull_polygon, axis=1
)

# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()

# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']

# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0]  # The first principal component represents the longer axis

# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis

# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj

# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3

# Classification function
def classify_region(projection, threshold_1, threshold_2):
    if projection <= threshold_1 or projection >= threshold_2:
        return 'pole'
    else:
        return 'center'

# Apply classification on filtered points
adata_filtered.obs['region'] = [
    classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')

# Visualize results
sq.pl.spatial_scatter(
    adata_filtered,
    size=10,
    library_id="spatial",
    figsize=(5, 5),
    frameon=False,
    shape=None,
    color=["region"]
)
adata_filtered.write('dyno4.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  _cax = scatter(
No description has been provided for this image
In [19]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA

adata1=ad.read('dyno_dm_5.h5ad')

# Extract spatial coordinates
points = adata1.obsm['spatial']

# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)

# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
    point = Point(row['x'], row['y'])
    return point.distance(polygon.boundary)

# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
    compute_distance_to_boundary, polygon=hull_polygon, axis=1
)

# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()

# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']

# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0]  # The first principal component represents the longer axis

# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis

# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj

# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3

# Classification function
def classify_region(projection, threshold_1, threshold_2):
    if projection <= threshold_1 or projection >= threshold_2:
        return 'pole'
    else:
        return 'center'

# Apply classification on filtered points
adata_filtered.obs['region'] = [
    classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')

# Visualize results
sq.pl.spatial_scatter(
    adata_filtered,
    size=10,
    library_id="spatial",
    figsize=(5, 5),
    frameon=False,
    shape=None,
    color=["region"]
)
adata_filtered.write('dyno5.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  _cax = scatter(
No description has been provided for this image
In [ ]:
 
In [20]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA

adata1=ad.read('dyno_dm_6.h5ad')

# Extract spatial coordinates
points = adata1.obsm['spatial']

# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)

# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
    point = Point(row['x'], row['y'])
    return point.distance(polygon.boundary)

# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
    compute_distance_to_boundary, polygon=hull_polygon, axis=1
)

# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()

# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']

# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0]  # The first principal component represents the longer axis

# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis

# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj

# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3

# Classification function
def classify_region(projection, threshold_1, threshold_2):
    if projection <= threshold_1 or projection >= threshold_2:
        return 'pole'
    else:
        return 'center'

# Apply classification on filtered points
adata_filtered.obs['region'] = [
    classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')

# Visualize results
sq.pl.spatial_scatter(
    adata_filtered,
    size=10,
    library_id="spatial",
    figsize=(5, 5),
    frameon=False,
    shape=None,
    color=["region"]
)
adata_filtered.write('dyno6.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
  _cax = scatter(
No description has been provided for this image
In [ ]: